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The capture of compact bodies by black holes in galactic nuclei is an important prospective 
source for low frequency gravitational wave detectors, such as the planned Laser Interferometer Space 
, Antenna. This paper calculates, using a semirelativistic approximation, the total energy and angular 

■ momentum lost to gravitational radiation by compact bodies on very high eccentricity orbits passing 

' close to a supermassive, nonspinning black hole; these quantities determine the characteristics of 

the orbital evolution necessary to estimate the capture rate. The semirelativistic approximation 
improves upon treatments which use orbits at Newtonian-order and quadrupolar radiation emission, 
and matches well onto accurate Teukolsky simulations for low eccentricity orbits. Formulae are 
presented for the semirelativistic energy and angular momentum fluxes as a function of general 
orbital parameters. 
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' Proposed space-based gravitational wave interferometers such as the Laser Interferometer Space Antenna (LISA) 
lO , will have good sensitivities in the low frequency gravitational wave band, from about 10~^ Hz up to about 1 Hz. A 
' promising source of gravitational waves in this band is the extreme-mass-ratio inspiral (EMRI) of compact objects 
, (stellar mass black holes, neutron stars, white dwarfs and even main sequence stars) into massive black holes. Current 
^ '■ estimates suggest we might detect as many as a thousand EMRI events over the course of the LISA mission 
^ \ Gravitational waves from extreme mass ratio captures will serve as a direct probe of the innermost population of 
bJO- compact objects around galactic central black holes, and also provide information on the growth history of such black 
holes out to a significant redshift (z 1). In addition to probing the stellar population of galactic nuclei, gravitational 
waves from EMRIs provide a map which encodes the geometry and structure of the black hole spacetime [2] , allowing 
a direct comparison of the astrophysical black hole to the black hole solutions of general relativity. This technique 
has been called "holiodesy", in analogy with satellite geodesy, which observes the motion of small satellites around 
the Earth to map out the structure of the planet's gravitational field. 

The computation of gravitational radiation from stellar orbits has a long history. The classic treatment was that 
of Peters and Mathews 0] , who computed the gravitational wave emission from stars on purely Keplerian orbits 
in flat space. The problem of generation and propagation of gravitational waves in a Kerr background was addressed 
by the work of Teukolsky and Press in the early 1970s [1, 0, who developed a perturbation formalism in the Kerr 
background. Subsequent work using this formalism has progressed to the point where the emission from a particle 
on any orbit in the Schwarzschild spacetime [8] or on a circular inclined Q or eccentric equatorial [To*! orbit in the 
Kerr background can be treated. Computing the inspiral of stars on eccentric non-equatorial orbits in Kerr required 
overcoming some technical obstacles [9[ , but first results are now available [ll|, ■ 

In this paper we compute a star's orbital trajectory by solving the geodesic equations of motion around the black 
hole, rather than using Keplerian orbits. Exact geodesies of the Schwarzschild spacetime are considered, particu- 
larly orbits with high eccentricity, including marginally bound (parabolic) and unbound (hyperbolic) orbits. Here 
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"parabolic" is a statement about the orbital energy E which labels the geodesic trajectory, rather than a statement 
about the geometric shape of the orbit in Euclidean geometry. As will be seen in later sections, the orbital trajectories 
around black holes can exhibit "zoom-whirl" behaviour, looping around the black hole more than once (i.e., the change 
in azimuthal angle Acj) > 27r) on any given orbital pass. Even in these situations one can still speak of quantities, 
such as the eccentricity of the orbit, which correspond in some formal sense to their Keplerian counterparts. Taking 
these geodesic trajectories as the orbits of the source, we approximate the gravitational radiation using the classic 
quadrupole formula at Newtonian-order. 

This method of obliging the orbiting body to follow a geodesic of the spacetime, while using the quadrupole 
approximation to calculate the gravitational wave emission, has been termed the "semirelativistic approximation" 
by its originators, Ruflini and Sasaki p^ . Experience with more accurate, although computationally intensive, 
perturbative calculations has shown that when particles make close approaches to the central black hole, in particular 
when they are relatively close to the unstable circular orbit (UCO) of the potential (see III B] for a definition), the 
properties of the relativistic gravitational potential are of critical importance in determining the gravitational wave 
flux. As will be seen, employing a more exact description of the particle trajectory (spacetime geodesies) together 
with an approximation of the wave flux is much more accurate in these cases than using a consistent Newtonian-order 
approximation. 

The semirelativistic approach complements the more complex Teukolsky based computations in several ways. First, 
technical difficulties and the demands of computing power have made Teukolsky calculations computationally difficult 
for orbits with eccentricity near unity, a regime where the work presented here is designed to work well. Second, because 
the computationally intensive Teukolsky approach is not practical for use in conjunction with typical simulations of 
the clusters of stars in galactic centers and their capture rates by the central black hole, it is useful to look for more 
convenient approximate methods, which are sufficiently accurate for reliable results. 

An extension of the semirelativistic approach is currentl y fi nding use in the computation of approximate EMRI 
waveforms for use in scoping out LISA data analysis Iisl. [iGj . A particle inspiral trajectory is computed by 
integrating post-Newtonian expressions for the energy and angular momentum fluxes. Integration of the Kerr geodesic 
equations along this trajectory yields the particle position as a function of time, from which a waveform is computed 
from an approximate quadrupole moment tensor generated as in the semirelativistic approach. These "numerical 
kludge" waveforms are inconsistent in that the energy and angular momentum content of the waveforms differs from 
the change in energy and angular momentum of the particle orbit which is nominally emitting the radiation. This 
energy inconsistency means that some of the results that have been obtained using kludged waveforms, such as signal- 
to-noise ratios [l|, are inaccurate. The results of the semirelativistic calculations presented here allow us to estimate 
the magnitude of this inconsistency for inspirals into Schwarzschild black holes (see section UlI C|) . In the future, using 
semirelativistic fluxes (once these are extended to spinning black holes) in the kludge in place of the post-Newtonian 
fluxes employed currently will yield consistent inspirals. 

The key results of this paper are summarised as follows: 

1. Numerical results are presented which enable us to explore and evaluate several approximate methods of cal- 
culating energy and angular momentum fluxes from EMRI orbits, and therefore the evolution of those orbits 
(see section ini|. Signiflcantly better evolutions (compared to the results of exact Teukolsky-based calculations) 
can be built out of the semirelativistic formalism developed here, but improved orbital evolutions can also be 
obtained from classic gravitational radiation estimates (Peters and Mathews) simply by choosing to work with 
"geodesic parameters" instead of "Keplerian parameters", with little consequence to computational cost. See 
section UlTSI and Figd 

2. Analytic expressions are derived for the energy {AE) and angular momentum (AL^) radiated in gravitational 
waves for a single orbital pass near a black hole, as a function of the orbital parameters, which exactly reproduce 
our numerical results. See section lTlI D 2l and Eqs. ([55|) and ([57)) . for the case of parabolic orbits, and appendix IA II 
for a more general discussion. 

3. Fitting functions are given which reproduce approximately, but to high accuracy, the analytic and numerical 
results. These functions are relatively simple expressions which could be conveniently used in place of consistently 
Newtonian-order expressions such as those of Peters and Mathews, which are significantly less accurate for orbits 
with very close periapses (see section [III D H and appendix IA 2[) . Although we present fits to the semirelativistic 
results only, the fitting functions have more general applicability and it should be possible to derive a fit of the 
same form to Teukolsky-based results once these are available for generic orbits. 

The remainder of the paper will be organised as follows. In section [TT] we describe the semirelativistic scheme which 
we use to model the gravitational radiation from EMRIs. In section|TTT]we present fluxes calculated using this approach 
and compare these with more accurate Teukolsky-based results, as well as with the consistently Newtonian results of 
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Peters and Mathews. We also present analytic formulae which reproduce our numerical results, and discuss the case 
of hyperbolic orbits which are initially unbound but become bound to the black hole via gravitational bremsstrahlung. 
Finally in section [IV] we summarise our most important results and findings. 

Throughout this paper, geometric units where G = c = I are employed unless otherwise specified. 



II. MODEL OF GRAVITATIONAL RADIATION 



A. Quadrupole approximation 

The energy and angular momentum carried away by gravitational waves from a weak-field, slow-motion source can 
be computed using the quadrupole formula p7j 

dE 1, , 

^ = -^{l,.I,.), (1) 

and 

= —-^^jkliikalal) , (2) 

where summations are implied over repeated indices, ejki is the permutation symbol, and Xjk is the reduced quadrupole 
moment of the system, 

Ijk= J p{xjXk - ^6jkr'^)d^x . (3) 

The angle brackets ( ) in ([l])-(l2|) indicate averaging over several orbits, but parabolic trajectories (our main focus) 
do not have periodic orbits. Indeed, the period of a parabolic orbit is infinite, so the average energy flux over the 
whole orbit is zero. Therefore it is convenient to work instead in terms of AE and AL^, the total energy and angular 
momentum radiated over a single orbit, which are in general finite. 

The corresponding gravitational waveform in transverse-traceless gauge is given by jl7j 



. TT _ 2 
jk — ~ 



Pjl -^Im Pmk 2 Pjf^ Pynl -^Irn 



Pjk = Sjk - nj Uk, (4) 



in which rij denotes the direction of propagation of the wave and r is the proper distance to the source. 

For orbits in the weak-field, far from the black hole, the quadrupole formula applies, the source particle orbit is 
Keplerian, and the radiation field reduces to the Peters and Mathews result. For orbits that pass close to the black 
hole, the particle's geodesic orbit is no longer Keplerian and the motion is neither weak-field nor slow motion, so the 
quadrupole formula does not describe the wave emission precisely. As described above, correcting the emission formula 
requires use of black hole perturbation theory (Teukolsky methods) , which is computationally very challenging, but 
a significantly improved approximation can be obtained by using the quadrupole formula with source particle orbits 
modified to be a geodesic of the black hole spacetime. 

To do this, first identify the Cartesian coordinates, Xj, in the quadrupole moment expression ^ with coordinates 
in the Schwarzschild spacetime. Treating the Schwarzschild coordinates {r, 6, 0} as spherical polar coordinates, define 
a set of pseudo-Cartesian coordinates by 

= (r sin6' cos r sin6' sin0, r cos6') . (5) 

In these coordinates, one can solve the Schwarzschild geodesic equations (see section IIIBI) to compute the particle 
orbit Xj(t), and substitute the resulting trajectory into the quadrupole moment expression ([3]) to compute 



jjk ^ 
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(6) 



where m is the mass of the particle. Finally, we compute an estimate of the energy and angular momentum radiated 
using expressions dJl-dl]). 

This approximation for gravitational wave emission was first applied by Ruffini and Sasaki, who termed it a 
"semirelativistic approximation" 13], since it makes use of the fully relativistic orbit, but only a weak-field expression 
for the gravitational waves. The approach is equivalent to attaching the compact body to a string in fiat space and 
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forcing it to move on a path that corresponds to a geodesic of the Schwarzschild potential. In reahty, the inspiralling 
body does not follow a geodesic, due to the effect of radiation reaction on the orbit. The loss of energy occurs 
continuously, so particle trajectories depart from a true geodesic path continuously. Instead of stable orbits, particles 
follow inspiralling paths with a steadily decreasing average radial distance from the center. However, in the typical 
case for extreme mass ratio inspirals, in which the rate of energy loss per orbit is small, the actual particle trajectory 
looks similar to a geodesic orbit for long periods, so one is justified in making an adiabatic approximation Q; simply 
assume the body evolves through a sequence of geodesies and determine this sequence using the energy and angular 
momentum fluxes from each geodesic orbit. 

The adiabatic approximation will break down when the orbital parameters change significantly on the timescale of 
a single orbit, which occurs only very close to the final plunge. The trajectory and waveform in this region must be 
computed using the computationally intensive self- force formalism (see [l^ for a review). 



B. Geodesies 

The equations governing geodesic motion in the Schwarzschild spacetime, in the usual Schwarzschild coordinates, 
are given by 
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^Hdfj^^- (8) 

where r is the proper time along the geodesic, is the conserved specific angular momentum of the particle, E is 
the conserved specific energy and M is the mass of the central black hole. We have taken advantage of spherical 
symmetry to assume an equatorial orbit, = 7r/2, without loss of generality. The radial equation of motion ([7]) may 
be written as a cubic polynomial divided by . The cubic has one, two or three roots depending on the values of 
E and L^- These roots correspond to turning points of the radial motion. Orbits with a single turning point plunge 
into the black hole and correspond to energies 

E^-l > ^ - — - "^V" , (10) 



where 





2M 


2MLl 


rj - 


ri 
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For bound orbits, it is possible to define an orbital eccentricity by analogy with the Keplerian case. Define the position 
of the apoapse of the orbit to be 

'^a = |^rp, (11) 

where is the radius of the per iapse. Equation pip is used to define the eccentricity of a geodesic in terms of the 
turning points of the orbit (sl. Il0j|. This definition carries over to parabolic {E^ = 1) and hyperbolic orbits (E"^ > !)• 
In the parabolic case, the radial geodesic equation ([7]) has only two turning points (the apoapse is 'at infinity'), but 
definition (fTTj) holds with e = 1. In the hyperbolic case, one of the turning points has r < 0; using this in (fTTj) one 
finds e > 1, and so in this case we call that turning point the apoapse. 

For this definition for eccentricity, we use the parameters (fp,e) to characterise the orbit, instead of (E,Lz). The 
energy and angular momentum are related to the periapse and eccentricity by 



p- H I M(l-e)(4Af-(l + e)r,) 

^ rp((l + e)rp-(3 + e2)M) ' ^ > 



(1 + e) Tp 



(l + e)^-(3 + e2) 



(13) 
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The radial geodesic equation becomes 



v^-l)(^-j {ra~r){r-rp){r~r^) (14) 

where the apoapse, r^, and energy E are given by equations (jlip and (|12p . and the third root of the potential is given 
by 

. - 2(1 + e)., ^^^^ 



((l + e)rp-4M) ' 

For any given eccentricity, there is a minimum value for the periapse below which the orbit plunges directly into the 
black hole. This occurs when the two inner turning points of the geodesic equation, r_ and r^, coincide. A geodesic 
with precisely this periapse asymptotically approaches a circular orbit as it nears the periapse, and spends an infinite 
amount of time whirling around the black hole. The asymptotic circular orbit is an unstable orbit of the gravitational 
potential, and we will refer to it as the 'unstable circular orbit' (UCO). The radius of the UCO determines the 
minimum periapse for geodesies of a fixed eccentricity. Equating r_ and Tp yields an expression for the UCO in terms 
of e 

rc/co^^^^M. (16) 
1 + e 

The statement that orbits with Tp < rjjco ai'e plunging is equivalent to the relationship (|10p between the energy and 
angular momentum (see |8] for an equivalent relation in terms of the semi-latus rectum) . If e = the UCO is at the 
familiar innermost stable circular orbit, r = 6M. In the extreme hyperbolic limit, e —> oo, the UCO approaches the 
horizon r = 2M . Parabolic orbits (e — 1) have a minimal periapse of r = 4M. Cutler, Kenncfick and Poisson [8] also 
discuss the UCO, but they call the line rp = rjjco the 'separatrix', since it separates bound from plunging orbits in 
phase space. 

These orbital properties can be understood by considering the radial gravitational potential V{r,Lz), which is 
illustrated in Figure [T] for a typical zoom- whirl orbit. The characteristic feature of these highly relativistic potentials 
is the maximum at rjjco- As an inspiral approaches plunge, the orbit is within the potential well but close to the top 
of the well. That is, the periapse lies close to the UCO. The particle thus zooms out to apoapse and back, but loiters 
close to periapse, whirling several times around the black hole on a nearly circular orbit before zooming out to apoapse 
again. As one approaches the UCO, the more exaggerated the whirl phase gets and the closer the resemblance to an 
unstable circular orbit. 



C. Waveform structure 



The waveforms resulting from zoom-whirl orbits are easy to comprehend. During the long apoapse passage the 
motion of the source is relatively slow, and the amplitude and frequency of the gravitational wave produced are both 
low. Near periapse the motion is much more rapid and the signal has much higher amplitude and frequency. If the 
whirling phase of the orbit consists of one or more complete revolutions about the central body then the waveform 
will have several cycles (two for each revolution) . The result is a waveform with a very low frequency (determined by 
the radial period of the orbital motion) and low amplitude superposed with a burst of short duration (relative to the 
overall period) and relatively high amplitude whose frequency is determined by the azimuthal period of the orbital 
motion. An example waveform is shown in Figure [21 corresponding to the orbit indicated in Figure [TJ 

Note that while the radial frequency is much too low for detection by LISA, the azimuthal (0) frequency does fall 
in the LISA bandwidth for the orbits of interest. Although there is a low probability of detecting these bursts since 
they are too brief and infrequent (typically occurring once every few years or even longer, depending on the radial 
period) to have high signal-to-noisc, the background of all such bursts occurring throughout our neighbourhood of 
the Universe will create an astrophysical background of noise from which other sources must be subtracted [19]. 



III. ENERGY AND ANGULAR MOMENTUM FLUXES 



The semirelativistic approximation is constructed by integrating approximate rates of energy and angular momen- 
tum flux over relativistically accurate geodesic orbits. A consistent approximation would require that the particle 
orbit be approximated to the same level of accuracy as employed for the fluxes. A Newtonian-order approximation. 
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FIG. 1: Radial gravitational potential for a zoom-whirl orbit. The dashed line corresponds to the energy of the orbit. The 
orbit oscillates in the region where the potential (solid curve) lies below the energy line. If the energy is too high and the orbit 
passes inside the maximum of the radial potential (ruco), the particle plunges into the black hole. 



such as that of Peters and Mathews, makes use of Keplerian elliptical orbits in flat spacetime and the fluxes are of 
quadrupole order only. There might not appear to be much justification for using accurate orbital paths but retaining 
approximate fluxes. For orbits which never come close to the central black hole the semirelativistic scheme does not 
improve signiflcantly on fully consistent Newtonian approximations - in nearly flat regions of spacetime all reasonable 
approximations fare well. However, orbits with small periapse distances are a very different case. More accurate 
schemes (such as those based on solution of the Teukolsky equation) show that the radiation from orbits which come 
close to the black hole show features that are greatly modified by the strongly curved spacetime and which are quali- 
tatively different from those seen at large radii from the black hole. Many of these features arise from the properties 
of the geodesies in the strong-field regime and therefore, as argued in such features can be modelled by schemes 
which combine exact geodesies with approximate fluxes. This approach shows significant improvements over the weak 
field approximation [3|, as we will see in the next section, while being considerably less expensive computationally 
than solving the Teukolsky equation. 



A. Comparison to Peters and Mathews and Teukolsky results 

The expressions derived in ^] for the energy and angular momentum fluxes from a Keplerian orbit are 

/d^\ ^ ^32^ (l-e)l / 73 37 \ /r^X- 

\dt/ 5 M2 24 ^ 96 j VmJ ^ 



dt / 5 M (i + ef 



These are equations (5.4) and (5.5) of 0, but written to lowest order in the mass ratio, m/M, in the extreme-mass- 
ratio limit M = mi ^ m2 = m. The eccentricity and periapse of a Keplerian orbit are given in terms of the energy 
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FIG. 2: Sample gravitational waveform (top) from the zoom-whirl orbit indicated in Figure [T] We show the plus polarisation 
of the gravitational wave as a function of time. The radiation is emitted predominantly in a high frequency burst during the 
whirl phase of the orbit. For comparison, the waveform from the Keplerian orbit with the same parameters is shown in the 
bottom diagram. 
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and angular momentum by 

To use pT)) - (IT8l) in the strong-field regime, the natural way to proceed is to evaluate the fluxes in equations (|17 p - (fT8|) 
for the Kcplcrian orbit with the corresponding energy and angular momentum, i.e., substitute and from (jl9p 
into (fT7| -(fT8 |) ("Peters and Mathews with Keplerian parameters"). This approach runs into difflculties however, since 
Keplerian orbits do not exist for certain valid choices of E and L^, for example if > A/^/(l — E"^) the Keplerian 
eccentricity is undefined. An alternative way to proceed is to compute the geodesic eccentricity and pcriapse using 
expressions ([T2|) -([T3 l) and use these in the flux formulae p7 |) -(fT8 |) . thus identifying geometrically similar orbits ("Peters 
and Mathews with geodesic parameters"). 

In Figure [3] we compare the fluxes computed in these three ways: Peters and Mathews fluxes using Keplerian 
parameters, Peters and Mathews using geodesic parameters, and the semirelativistic approximation, all as a function 
of geodesic (relativistic) periapse for flxed geodesic (relativistic) eccentricity of e = 0.99. For large periapse, the three 
approximations agree as expected, but once the periapse becomes moderate {vp < 5QM), the Peters and Mathews 
expression with Keplerian parameters begins to differ quite significantly from the other approximations. In the strong- 
field region (rp < lOM), the semirelativistic approximation begins to differ significantly from both applications of the 
Peters and Mathews formula, predicting greater fluxes of both energy and angular momentum. 

To verify that the semirelativistic results are an improvement over Peters and Mathews, rather than merely being 
different, the approximation can be compared to perturbative results from integration of the Teukolsky equation. Very 
few results are available for high eccentricities in the Teukolsky formalism, so the comparisons here are shown at lower 
eccentricities. In Figure IH the semirelativistic and Peters and Mathews fluxes are compared to Teukolsky calculations 
[ifll] for orbits with eccentricity e = 0.5 and a variety of periapses. As one would expect, the semirelativistic approx- 
imation is not superior to a consistent Newtonian approach for periapses greater than about ^ 50M (sometimes it 
does worse and sometimes better than the Peters and Mathews results, but never extremely different). For periapses 
less than ~ 50 Af, the semirelativistic method improves significantly upon the consistent Peters and Mathews approx- 
imation. The improvement gained using the semirelativistic approximation was also noted by (20j in comparisons to 
a selection of Teukolsky-based results. Thus for the type of orbit of interest to this paper (highly eccentric orbits with 
close periapses), consistent Newtonian approximations should be regarded with suspicion, but approximations which 
make use of exact geodesies (like the semirelativistic approximation) , will fare very well qualitatively and quite well 
quantitatively, as long as the periapse is not extremely small. 

What is perhaps more surprising is that one may obtain an improved approximation from the Newtonian-order 
expressions (i.e., Peters and Mathews) if one carefully chooses the Newtonian parameters which are to be equated 
with the "true" curved space parameters (i.e., geodesic parameters rather than Keplerian parameters). While the 
semirelativistic approximation is always an improvement over this in the strong-field regime, the gain is only significant 
for very close periapses, Vp < lOM. In fact, for small eccentricities there is no significant gain using the semirelativistic 
fluxes. For a circular orbit of radius tq, the quadrupole formulae dl])-© tell us that {dE/dt) = —32rQQ^ and 
(dLz/dt) = — 32rgn^, where fi^ — d(j)/dt is the angular velocity. For both a Keplerian orbit and a circular geodesic 

3 

of the Schwarzschild metric, fl^ = l/r^ . Therefore, the standard Peters and Mathews result is recovered exactly for 
circular orbits using either geodesic or Keplerian parameters. 

We are primarily interested in highly eccentric orbits, for which the semirelativistic results are a significant improve- 
ment over any method based on Peters and Mathews. Nonetheless, if one does not wish the additional computational 
burden of evaluating more accurate semirelativistic flux expressions, a significant improvement can still be gained by 
evaluating the Newtonian fluxes using geodesic parameters. 

B. Phase space structure 

An inspiral sequence may be constructed from the semirelativistic fluxes by integrating {dE/dt, dL^/dt). While 
the duration of the inspiral depends on the value of dE/dt, the sequence of geodesies that the inspiral passes through 
in the (E,Lz) phase space depends only on the ratio dE/dLz- Equivalcntly, in the {rp,e) phase space, it depends 
only on 



9 




5 10 15 20 25 




FIG. 3: Comparison between the semirelativistic results and Peters and Mathews as a function of periapse for orbits with 
fixed geodesic eccentricity e — 0.99. The solid lines are the results from the semirelativistic approximation discussed here. The 
dashed and dotted lines are the Peters and Mathews results with geodesic parameters and Keplerian parameters respectively. 
We use a logarithmic vertical scale and plot the absolute value of the energy (upper plot) and angular momentum (lower plot) 
fluxes. 

It turns out that the semirelativistic approximation reproduces the ratio dE/dL^ to a very high accuracy when 
compared to the Teukolsky results. While the value of AE can differ by as much as 25%, the ratio dE/dL^ is 
recovered to better than 5%. This means that the structure of the semirelativistic phase space will be a good 
approximation to the true structure, although there is some error in the estimated duration of inspirals. 

This is an interesting result from the point of view of detection of EMRIs with LISA. An error in the timescale 
of an inspiral can be partially corrected by a change in the mass ratio. Since the phase space trajectory is well 
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FIG. 4: Comparison between accurate Teukolsky results and various approximations, for orbits with fixed eccentricity e — 0.5 
and a variety of periapses. The plots show the ratio of the flux computed using a given approximation to the flux obtained 
from the Teukolsky calculation. The upper plot shows the ratio of the energy fluxes, while the lower plot is the ratio of the 
angular momentum fluxes. In both plots, the solid lines are the semirelativistic results. The dashed and dotted lines are for 
Peters and Mathews, evaluated with geodesic parameters and Keplerian parameters, respectively. The latter lines cut off at 
small periapse since there are no corresponding Keplerian orbits in that region. 



approximated, an inspiral waveform computed using this approach might be a moderately good fit to a true inspiral 
waveform with a slightly different mass ratio, and therefore could be used as a detection template over sufficiently 
short time stretches. This may not be practical, since the error in dE/dt is not a constant factor, which would require 
varying the mass ratio over the inspiral. Moreover, other features that these approximations do not include (such as 
the "conservative" part of the self-force) may lead to rapid de-phasing of the kludge templates. This is nonetheless 
an interesting result. 

The accuracy with which the phase space structure is reproduced can be understood by considering what happens 
at the extremes of an inspiral. When the periapse is very large, the semirelativistic approximation is good, and is 
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expected to reproduce dE/dLz accurately. For an orbit that sits on the separatrix between plunging and non-plunging 
orbits (rp ~ rijco)j the geodesic spends an infinite amount of time whirling about the black hole on a nearly circular 
orbit at the UCO. The flux of energy and angular momentum is totally dominated by the circular part of the orbit. 
For any radiation field in a spherically symmetric spacetime, the energy and angular momentum fluxes carried away 
from a circular orbit obey the relation dE/dt = il.^{dLz/dt), where il^ is the angular velocity on that orbit [8.]. The 
quadrupole formulae (IT])-((2]) reproduce this result for a circular orbit. Since we use the correct input geodesic, the 
semirelativistic approximation should and does give dE/dLz accurately on the separatrix. Since we have the correct 
result in both extremes, it is perhaps less surprising that we also do quite well at points in between. 

Figure[5]illustrates some inspiral trajectories in the (rp, e) plane. The trajectory properties are determined largely by 
two curves - the separatrix where the inspiral ends as the particle plunges into the black hole and the locus de/dt = 0. 
These are both marked on Figure [S] In general an inspiral will begin with high eccentricity and moderate periapse. 
Both the periapse and eccentricity initially decrease, and this evolution continues until the trajectory intersects the 
de/dt = curve. After that point, the periapse continues to decrease, but the eccentricity increases until the trajectory 
reaches the separatrix and the particle plunges. As expected from previous arguments, these general properties are in 
good agreement with results based on Teukolsky calculations [s"] and quite different to Peters and Mathews inspirals 
(which, for instance, have monotonically decreasing eccentricity). The location of the separatrix is a property of the 
geodesies, and is therefore precisely reproduced in this approximation. The de/dt — curve depends on the expression 
used for the energy and angular momentum fluxes and is different here, but only slightly. In this approximation, the 
de/dt = curve intersects the e = axis at rp — 6.770 M, compared to Vp — 6.681 M in the Teukolsky case 8]. 

The increase in eccentricity prior to plunge is a generic feature of EMRIs, but it is as much a property of the 
radial potential as it is of the flux model. As discussed earlier, realistic gravitational waves will give E — fl^ — 
((1 + e)/(2 (3 + e)))i Lz on the separatrix. The leading order piece in both the numerator and denominator of 
equation (|^D|) thus vanishes in the limit rp — > rucois), but the leading correction to both is from the logarithmic 
piece of dE/dLz, and hence we find drp/de < 0. However, this conclusion still holds if the fluxes do not satisfy the 
circular orbit condition and the cancellations do not occur. The coordinate derivatives are such that, independent of 
the value of dE/dLz on the separatrix, 

drp ^ 4(3 + e) 

de ^ (l-e)(l + e)2 ^ ^ 

The nature of the potential thus forces either rp or e to increase in the approach to plunge. 

A final point to note from Figure [5] is that e oc e near e = 0. This property of the inspirals means that an initially 
eccentric orbit cannot circularise in this model, although the eccentricity at plunge can be arbitrarily small. The 
property once again derives from the circular condition E = Qj, Lz, which ensures that circular orbits remain circular 
under radiation reaction. This is discussed in more detail in [l5| , where corrections are given to enforce this relation 
in kludged inspirals. In the semirelativistic waveform model, the condition is automatically satisfied and no correction 
is required. 

C. "Kludge" waveform inconsistency 

As mentioned in the introduction, waveforms based on the semirelativistic approximation are being used exten- 
sively to scope out LISA data analysis [l3, [3, [3] ■ The waveforms are generated by first constructing an inspiral 
trajectory and then using the semirelativistic construction of the quadrupole moment tensor to compute a waveform. 
In the most basic form of this "kludge" [l3|, the phase space trajectory for inspirals into non-spinning black holes 
is computed by integrating the Peters and Mathews energy and angular momentum fluxes P7)l . (fT^ . This leads to 
an inconsistency since the energy and angular momentum content of the gravitational waves differs from the energy 
and angular momentum lost by the inspiralling particle that is nominally generating the radiation. We can estimate 
this inconsistency using the semirelativistic results. A phase space trajectory is generated using the fluxes ^T7\ . 
and equation (^0]) . We choose to specify the eccentricity of the inspiral at plunge and integrate backwards along the 
inspiral trajectory. By integrating the semirelativistic fluxes along this trajectory, we calculate the total energy and 
angular momentum content of the gravitational waves. Figure |6] shows the ratio of the gravitational wave energy flux 
to the change in energy of the particle orbit as a function of the time until plunge (in units of M'^/m). Time along 
the inspiral trajectory therefore decreases toward the right. There is a curve for each eccentricity at plunge from 
e = 0.1 to e = 0.9 in intervals of 0.1. We see that there is a significant inconsistency in the kludge waveforms. For 
low eccentricity at plunge, the kludge gravitational waves contain less energy than they should, but for eccentricity 
at plunge greater than about 0.25, they contain too much energy, as much as a factor of three in extreme cases. This 
means that signal-to-noise ratios (SNRs) computed from these waveforms are likely to be overestimates of the true 
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FIG. 5: Sample inspiral trajectories in the (rp,e) plane. The solid lines illustrate inspiral trajectories. The dotted line marks 
the separatrix and all points to the left of this line are plunging orbits. The dashed line is the locus of points with de/dt = 0. 
In the region between this line and the separatrix, de/dt > on all trajectories. 

SNRs. It is not clear from these results whether this discrepancy will be larger or smaller when the central black 
hole is spinning, but this will be investigated in the future [l3|. However, it is important to be aware of the existence 
and magnitude of this problem when interpreting results based on the kludge waveforms. If semirelativistic fluxes 
were used to integrate the phase space trajectories, there would be no such inconsistency and this might therefore be 
another future application of these results, once they are extended to spinning black holes. 

D. Analytic results 

The previous results have shown the usefulness of the semirelativistic approximation, but the described method, 
based on integration of the geodesic equations, is not easy to implement in numerical simulations. In this section, 
we present some analytic results based on the semirelativistic approximation which can be easily evaluated without 
numerical integration of orbits. 

1. Fitting functions for AE and ALz 

A useful tool for simulations is a fitting function that has a simple form and which reproduces the semirelativistic 
results for AE and AL^ to reasonable accuracy. For a geodesic of given eccentricity, the periapse can have any value 
between the UCO for that eccentricity ([TB)) and infinity. For large values of the periapse, the orbit is entirely within 
the weak-field region of the spacetime. The orbit and radiation will therefore look very like those from a Keplerian 
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FIG. 6: Ratio of the energy content of "kludge" gravitational waves to the change in energy of the source particle orbit, relative 
to the energy at plunge. This is shown as a function of time until plunge (in unite of Af^/m) for eccentricities at plunge 
Epi — 0.1, 0.2. ..0.9 (from lowermost curve to uppermost). 



orbit, as described by expressions p?)) - p^ Multiplying these expressions by the Keplerian orbital period, the 
energy and specific angular momentum lost on a single pass may be found to be 
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AL, = -51^,,.,^ + (^) . (23) 



In the parabolic case (e = 1), these become 



AL, = — Gttto 



As the periapse approaches the UCO, the energy and angular momentum lost per pass increases. In fact, ignoring 
radiation reaction, the energy and angular momentum losses diverge as rp rjjco- This is because the geodesic with 
''p = ^uco spends an infinite amount of time whirling around the black hole with r w rjjco- In practice, radiation 
reaction will prevent this situation arising (for a discussion of the transition from inspiral to plunging orbits see [2l|). 
However, the energy and angular momentum lost should increase rapidly as the periapse approaches rucOi since the 
orbit whirls around the black hole an increasing number of times. 

As discussed previously, the whirl behaviour is a property of the geodesies. Thus, although the semirelativistic 
treatment of the radiation is only approximate, one will still see this divergent behaviour as rp — > r[/co, since the 
source for the radiation is an exact geodesic trajectory. This feature is missing in the Keplerian treatment [3[. For an 
extreme zoom and whirl orbit, most of the gravitational radiation is emitted during the whirl phase, when the particle 
is on an approximately circular orbit. It is reasonable to guess that the total losses due to gravitational radiation are 
roughly proportional to the number of times the particle whirls around the black hole. 
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In the parabolic case, one can estimate the number of whirls by computing the proper time taken for the orbit to 
pass from periapse to some "whirling radius" , r„ . This is found using to be 
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(26) 



3 /2 

The radius does not change significantly during the whirl phase, so approximating the numerator by 

' .|cosh-^2(^^-^^) 



V2M 



rp{rp-AM) 



2{rp -2M) 



(27) 



Using the same assumption, d^/dr w L/r^ and we estimate that while r < Vyj, the number of azimuthal cycles that 
the particle completes is 
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The radius should be chosen to define the start and end of the whirl phase. Our objective is to guess a functional 
form that approximates the energy and angular momentum loss when w tuco ^nd we assume that dE and dL^ are 
proportional to (|28p in that limit. This is likely to be a particularly good approximation for highly eccentric orbits, in 
which the 'zoom' and 'whirl' phases are quite distinct. Appropriate fitting functions should approach ((22l) and ((23l) in 
the limit rp — > oo and should diverge like (|28p as rp — > rjjco- Working once again in the parabolic case, the simplest 
such function is 
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(29) 



One could fix all three coefficients by matching the behaviour in the limit rp rjjco j but (j29p will not then necessarily 
reproduce the asymptotic form ([M)) . Instead, we fix and using an expansion near rp — ruco and then fix 
= -(857r/(12V2) + 64:A^V2B^) to match ^ asymptotically. 
The next section will demonstrate how exact expressions for the energy and angular momentum radiated in our 
model may be obtained in terms of elliptic integrals. Using these full analytic expressions, we can predict the values 
of the fitting function coefficients 
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= 0.752091, 



The equivalent fitting function for the angular momentum lost is 
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with coefficients 

A-^ = -i^ = 
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C^^ = - 67r 



A^' V2B^] = -4.149103 . 



(30) 



(31) 



(32) 



The fitting function (j29p can be used to match the lowest order terms in an expansion of AE near rp = ruco and 
Tp — > cx). It is possible to add additional terms to give a more general function which can match AE at arbitrary 
orders 
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In this, we fix Ne = 7 to give the correct leading order behaviour ([24)1 as rp — > oo. The parameter N indicates 
the order of the fit, i.e., the number of terms we include. The second and third series have terms in common, but 
writing the expansion in this way allows one to read off consecutive coefficients easily. The next section will show that 
expanding AE about the separatrix gives terms in (r^ — AAiy ln(rp — 4Af ) and in (rp — 4M)'^, while an expansion as 

Tp — > oo gives terms of the form . The coefficient of the j — term gives Aq , then the fc = term gives 

and the / = term gives Cq . Continuing in this way, the j,k,l — n terms determine A^, and respectively. 
Thus, an expansion to order N will match the lowest + 1 terms in j, k and I. A similar fitting form may be used 
for ALz/M , but with Ne replaced by iV^^ = 4 (once again, to reproduce the correct leading order behaviour (|25p as 
Tp oo). Figure [7] illustrates how the fitting functions converge as the order of the fit, N , is increased. We see that 
as N increases, the fit improves at large radii, but initially gets slightly worse at moderate radii before converging 
there also. The iV = 2 fit is accurate to about one percent everywhere, so we include these parameters here 

A^ = -0.141421, Af = 0., -1.20797, = 0.752091, Sf = -103.215, 

= 727.515, = -4.63464, Cf = 69.1683, = -439.378 

Ag-- = -1.13137, A[' = 0., A^' ^ 0., Bq' = 1.31899, = -53.4491, 

b!^' ^ 29.7857, C^' = -4.1491, C^' = 25.4129, C^^ = 15.1726 (34) 

This fitting function was derived using simple arguments about how the energy and angular momentum lost behave. 
These arguments are valid in general for radiation that is produced by a body orbiting in the Schwarzschild potential 
and will apply to fluxes computed using the Teukolsky formalism. In a separate paper [2^ we derive a fit of this form 
to Teukolsky data computed for parabolic orbits in [23], which even for iV = 2 is accurate to < 0.2% everywhere. 

This simple fitting function is clearly a useful and accurate way to evolve EMRI orbits. In the case of arbitrary 
eccentricity, a similar type of fitting function can be derived, but the coefficients A^ etc. are now functions of 
eccentricity. In the semirelativistic approximation, the functions can be evaluated explicitly. This is discussed in more 
detail in appendix I A 21 It is more complicated to compute a fit to Teukolsky-based fluxes, since the coefficients in the 
expansion must be further expanded as functions of eccentricity. However, it should be possible to derive a reasonable 
fit using a polynomial ansatz, of the form suggested by the semirelativistic results. Once sufficient Teukolsky-based 
data is available, this fitting procedure will allow us to generate a comparatively simple analytic expression for use in 
computation of EMRIs. 

2. Exact expression 

As mentioned in the preceding section, it is possible to derive exact analytic expressions for the radiation loss 
predicted by our quadrupole approximation. This is possible because in any axisymmetric spacetime, the rate of 
energy loss at a given moment in time cannot depend on the absolute value of the (j) coordinate of the particle, since a 
shift 4> ^ (/) + 00 will leave the spacetime unchanged (note that the energy flux in a given direction will be dependent 
on the relative difference in </> between the source and the observer). In the quadrupole approximation used here 
to compute the gravitational radiation, dE/dt is given by the square of the third time derivative of the quadrupole 
moment tensor. This is a function only of the particle coordinates r{t) and (pit). By the axisymmetry argument, the 
expression for dE/dt can depend only on r, f, f, r', </>, cj) and (f>. For geodesies of the Schwarzschild spacetime (also 
for equatorial orbits in the Kerr spacetime), d^/dr and dt/dr are rational functions of r only and (dr/dr)^ = V{r) 
is a cubic or quartic polynomial in r. Any derived expression, such as dE/dr — (dE / dt) / (dr / dt) , will therefore be 
a rational function of polynomials in r and \/V{r). It is known [l^] that the integral of any rational function of 
polynomials in x and y, where is a cubic or quartic polynomial in x, can be expressed in terms of elliptic integrals. 
One can therefore write the energy and angular momentum radiated in closed form in terms of elliptic integrals. 

By substitution of the geodesic equations ([8]), ([9|) and (fT4| into (Il])-([3]), we may write dE/dt and dL^/dt as functions 
of r and then integrate over one orbit. In the parabolic case, the energy loss is found to be 
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h{v) = - 2?/ (27850061568- 83550184704 y+ 117662445984 y^- 102686941680 
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FIG. 7: Error using fitting functions to approximate the analytic expressions for the energy (upper plot) and angular momentum 
loss (lower plot). In each plot, the absolute percentage error in the fit is shown for fitting functions of various orders, = 
2,--. ,6. 
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and 

72(2/) = (-72901570560+ 274404834816 y- 424693524096 y2 
+378109481088 - 249480499840 + 154011967968 
-84437171728 y^ + 31689370996 y^ - 6231594434 / + 321950817 
+27462280 + 4073612 y" + 696320 y^^ + 40960 y") . 

In this, K and E are the complete Ehiptic integrals of the first and second kinds respectively, defined by [24l |2 



K(fc)= / ^ E(A:)=/ \/l-k^ sm^(f>d^. (36) 

V 1 - fc^ sin^ 



The corresponding result for the angular momentum lost is 
AL;, 64 



24249225 Tp- (r^ - 2A/)3 



(37) 



where 



and 



(y) = y (181817664 - 363635328 y + 245236248 y^ - 49673460 y^ 

-7833906 y^ + 2016105 y^ + 283252 y^ + 35896 y^ + 4120 y^) 



52 (y) = (71285760 - 324389184 y - 468548880 y^ - 277856496 y^ + 54521424 y"* 
+6181872 y^ - 1630457 y^ - 238086 y^ - 31776 y^ - 4120 y^) 

These exact expressions can be used to derive the fitting function described in the previous section. As Tp — s- oo, 
the argument of the elliptic integrals, y/2 M/{rp — 2 M) —^0. In a series expansion of the integrals about fc = the 
lowest five orders in k cancel and one successfully recovers 

As r — > rjjco — 4 M, the argument of the elliptic integrals a/2 M/ {vp — 2 M) — > 1 and the elliptic integrals diverge. 
Using [i^l and some algebraic manipulation, the asymptotic forms of the elliptic integrals as /c — > 1 are found to be 

K(fc) = In (1 - fc2) + 2 ln2 - i (1 - k^) In (l - fc^) + O (l - k^) (38) 
2 8 

E(fc) = 1 - i (l - /c^) In (1 - fc2) + J^ln2 - ^ (l - k^) +0 ((1 - k^ f ln{l - k^)) . (39) 

The asymptotic form of ((35|) -([37 l) as rp ^ 4 M is 

dX 
V 

In this, X refers to either E or L^/M . The values of the constants are 



i^.„,„(i-4)+„+0(i-4) ,40) 



1 „ / 16370483137 ln(2 \ , 

PE = — F, qE = 2 [ ^ , 41 

5V2 V 53542288800 V5 2%/2/ 

4%/2 „ /l613849^/2 „ ^ , , ,\ 

P.. = — ,..^2(^^^^-2V2 1n(2)j. (42) 



^ AX = -A;, In - 4) + In (26;,) + O - 4) (43) 



The fitting functions (P^ - (PT|) may be similarly expanded near Vp = rjjco = 4 M 

M 
m 

In this, X once again refers to either E or Lz/M, and A''^; = 7, A^l^ = 4. Equating ([^0]) and one obtains the 

coefficients of the fitting function given earher ([50]) . ([5^ . 

A similar analysis can be performed for orbits of arbitrary eccentricity, and is described in appendix I A II The exact 
expressions are somewhat cumbersome and we recommend using the fitting function in most applications, since this 
performs extremely well. The exact expressions have been included for completeness, and to help explain why the 
fitting function works. 
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E. Hyperbolic captures 

In this paper we mostly focus on parabolic orbits, which serve as a useful model for all orbits which are likely to lead 
to sources of interest to LISA, since such orbits will always initially have eccentricities very close to one. However, 
objects can also be captured from orbits with e > 1. In such cases the orbit is unbound, but may ultimately inspiral 
if it makes a close approach to the central black hole and loses enough energy and angular momentum in doing so to 
become bound. Our results suggest that if the angular momentum of this orbit is low (close to the minimum =4), 
then the scattered body will become bound if it is on an orbit whose energy E is such that E"^ — 1 < m/M roughly 
speaking, where m/M ^ the mass ratio, is small. For larger angular momenta the amount of excess energy which can be 
radiated away on the first pass is smaller, and so the orbital energy must be even closer to unity for the body to become 
bound. Figure [8] shows which hyperbolic orbits can lead to captures, for low orbital angular momenta. The energy 
and angular momentum lost to gravitational waves by a particle on a hyperbolic orbit are given by the same equations 
(IA3p - (|A4p that apply to bound orbits, just by inserting e > 1 consistent with the definitions (fT^ - P^ ^. Figure [5] 
was generated by using equations p2| -(fT3 |l in conjunction with equation (|A3|1 to write /S.E = {m/M)FE{E,Lz) for 
hyperbolic orbits. Points on the curves obey the equation 

in 

E-1 = -—Fe{E,L,). (44) 

Fixing the energy, the angular momentum solution to (|44p is obtained by iteration. 

This figure indicates that there is another type of orbit which becomes bound after the first pass - those that are 
close to the separatrix. If the orbit is very close to the plunge line, it will also lose enough energy to become bound 
even if it has much more energy than E^ — 1 = m/M. The reason seems obvious after a glance at Figure [TJ If the 
energy of the scattered body is sufficiently close to 1 then it is close enough to the potential well in which bound 
orbits exist to lose sufficient energy on one pass to fall into the well. If the energy and angular momentum of the 
orbit are such that the particle's energy places it at the maximum of the potential, then the particle "whirls" around 
the central black hole at the radius of the potential maximum. The scattered body thus spends an abnormally long 
time near periastron and hence radiates an unusually large amount of energy, enough to become bound. Of course 
these bodies will generally plunge rather soon after capture because the amount of angular momentum radiated will 
also decrease the height of the potential barrier. In fact many of these orbits will plunge on their first pass, having 
dissipated enough angular momentum to shrink the potential barrier so they pass over it into the plunging region 
beyond. One has to keep in mind that our adiabatic approximation breaks down as one approaches the line separating 
stable from plunging orbits (depending on the mass ratio), so that we can give no definitive picture of what occurs 
in this regime except to say that the general behaviour is probably correct. Readers interested in the transition from 
inspiral to plunge should consult plj . 

Orbits that are scattered close to the separatrix line are "captured" , in the sense that they either plunge or become 
bound. Particles which are close to being parabolic orbits are also captured and may serve to modestly increase the 
capture rate for LISA. Particles passing very near to the black hole must thus pass between Scylla and Charybdis ^. 
If they pass too close to Scylla, through having energy only marginally greater than 1, then they are plucked from 
their unbound orbit by gravitational radiation reaction and end up in a bound orbit. If they approach with too great 
an energy for their angular momentum then they are sucked down by Charybdis and plunge into the black hole itself. 

IV. RESULTS AND DISCUSSION 

As with earth-based gravitational wave detectors like the Laser Interferometer Ground Observatory (LIGO), theo- 
retical predictions of source event rates and signal characteristics for LISA will play an important role in the successful 
operation of the observatory. Up to now, Newtonian-order estimates (like Peters and Mathews) have been widely 
relied upon to estimate waveforms and fluxes from extreme-mass ratio inspirals, even though much of what is of 
interest to LISA, even in the early stages of inspiral, occurs inside the region of relatively strong curvature close to 



^ This statement, as appeared in the published version of this paper, is not correct. The published erratum has been appended in 
Appendix IBI 

In The Odyssey Scylla and Charybdis were monsters who guarded opposite sides of a narrow straight through which ships must pass on 
their way from Troy in Asia Minor to Ithaca in Greece. Scylla was a six-headed reptile who plucked sailors from the decks of their ships 
as they passed, while Charybdis was an undersea monster which created a great whirlpool by its sucking in of seawater, the whirlpool 
causing ships to be pulled under if they strayed too near. 
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FIG. 8: "Hyperbolic" orbits (i.e., orbits with E > 1) tliat are captured after one close encounter with the black hole, for various 
mass ratios. Orbits whose energy and angular momentum place them above and to the right of the line for a given mass ratio 
remain unbound and are not captured. The line which begins at bottom left and curves around to top right in the figure is the 
separatrix line separating unstable plunging orbits (to its left) from stable orbits. 



the central black hole. The principle reason for this is simply ease of use. Even when accurate methods, such as the 
Teukolsky formalism and self-force calculations, prove capable of dealing with arbitrary orbits they may still be slow 
and cumbersome for many applications. This paper attempts to make available a range of techniques which combine 
ease of use with fairly robust accuracy over almost the whole inspiral of an extreme-mass-ratio binary. These results 
are of particular use for highly eccentric orbits, where fre quen cy domain Teukolsky calculations perform poorly (lo| 
and time domain codes have not yet been fully developed 23j . 

The key elements to take away from this study of the semirelativistic approximation are: 

• Simple analytic expressions to estimate the fluxes AE and AL^, suitable for use in computational endeavours. 

• The optimal choice of parameters with which to describe orbits which stray near the central black hole are 
the geodesic parameters {rp,e} rather than {E,Lz}'. the waveforms for orbits which have similar {rp,e} more 
closely match than orbits which have similar {i?, Lz} values. 



This second point cannot be stressed enough, as it applies to treatments which use the semirelativistic approximation 
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or Newtonian results; all approaches appear to be most accurate when the orbits are defined by the periapse distance 
Tp and eccentricity e rather than by energy and angular momentum. The reason is that when working in flat space 
relating the orbit to the curved-spacetime orbit with the same and e gives much better agreement with the curved- 
spacetime fluxes derived by exact methods (Teukolsky methods) than with the fluxes from the curved-spacetime orbit 
with the same E and Lz as the flat space orbit. This is one substantial improvement in accuracy (see Figure[3]) which 
can be made for no computational cost whatsoever. 

To gain further improvements, the fltting function (jAll[) described in section Fill D II can be used for only a small 
additional computational cost. Using the coefficients presented here, we have seen that it can accurately reproduce 
the energy and angular momentum fluxes computed using the semirelativistic approximation. However, it also has 
more general applicability. Once sufficient data has been obtained by numerical solution of the Teukolsky equation, it 
should be possible to derive a good fit to that data using the same fitting ansatz. This will provide a more practical 
expression for use in astrophysical calculations. 

In section IIII CI we made use of the semirelativistic results to estimate the inconsistency in "kludge" gravitational 
waveforms that are being used to scope out LISA data analysis [13, d, EB]- These waveforms are constructed in 
a similar way to the semirelativistic fiuxes described here, but the inspiral trajectory of the particle is computed 
independently of the waveforms using post-Newtonian results. We saw that the energy content of the gravitational 
waves can be as much as a factor of three greater than the energy lost by the particle orbit. This is an important 
point to bear in mind when interpreting results computed using these approximations. 

The semirelativistic formalism presented here should find uses in computational problems where speed is of concern 
(e.g., large numerical simulations) and the role played by the central black hole is important to the dynamics of 
individual particles in the problem. Such problems of interest might include new simulations of star cluster evolution 
in galactic nuclei to estimate the LISA EMRI event rate, or supermassive black hole inspiral simulations which seek to 
use interactions with stellar populations as a source of dynamical friction to bring the large black holes into proximity. 
In a companion pap er 22l| . we use the insight gained here, in conjunction with numerical results from solution of the 
Teukolsky equation [23|, to compute improved expressions for the inspiral timescale of capture orbits. The resulting 
expressions can be easily included in simulations of stellar clusters |27j, [2^ to improve estimates of capture rates. 

The semirelativistic approximation can also be applied to estimate energy and angular momentum fluxes from 
objects orbiting Kerr black holes. The procedure is more complicated due to the inclusion of spin and lack of 
spherical symmetry. In particular, it is not clear how to evolve the third integral of the motion, the Carter constant, 
for Kerr inspirals. However, by identifying Boyer-Lindquist coordinates with flat-space spherical polar coordinates 
and constructing the corresponding flat-space quadrupole moment tensor in the manner employed here, estimates for 
the energy and angular momentum fluxes from Kerr orbits may still be obtained. Preliminary results suggest that 
such semirelativistic estimates improve over standard post-Newtonian results [2^ for spinning black holes as well. To 
construct inspirals, the angular momentum and energy fluxes can be combined with "kludge" approximations for the 
evolution of the Carter constant This extension to Kerr will be described in a future paper. 
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APPENDIX A: ANALYTIC RESULTS FOR ARBITRARY ECCENTRICITY 

1. Exact expressions 

For orbits of arbitrary eccentricity it is also possible to derive exact expressions for the loss in energy and angular 
momentum, which reduce to psp and (|37p in the parabolic case. This is accomplished by writing the energy and 
angular momentum lost as a sum of integrals of the form 

/ (Al) 



By considering the derivative of \/ {ra — r){r — rp){r — r^)r/r" and using results in |2a l, we deduce 

2n - 3) A-P{1 - e)((l + e)rp - AM) / n-2\ 



/ n - 1 \ Af((l + e)rp - (3 + e^)M) {2n 

= 7; 7 ^'1-1 





(l-e)((l + e)rp-4Af) / 4eA/ 



(l + e)((l + e)rp-2(3-e)Af) \ V ((1 + e)rp - 2(3 - e)A/) 



M{{l + e)rp-4M) / (1 - e)((l -f e)rp - 4Af) 



(l + e)r3 y (l + e)((l + e)rp-2(3-e)Af) 




4eAf rp((l + e)rp-2(3-e)A/) ^ ( j 4e_ 



((1 + e)rp - 2(3 - e)A/) / (1 -f e)rp - 4Af \ V ((1 + e)rp - 2(3 - e)Af) 



.(A2) 



The functions K(A;) and E(fc) are the complete elliptic integrals of the first and second kinds ([55)1 . Using the recurrence 
relation (jA2p we can express the energy and angular momentum lost in terms of these elliptic integrals. We find the 



expression for the energy loss to be 




1673196525r6(l + e)^ ((rp 



16M 



^ 2M)((1 + (),■„- 2(1 - e)M))i 



X 



+ 



^(l + e)^-2(3-e) 



y(l + e)g-2(3-e)E 



+ < 



K 



( 




((1 + e)rp - 2(3 - e)M) 



((1 + e)rp - 2(3 - e)M) 




where 



/i(y,e)= 4608 (1 -e) (1 + e)^ (3 + e^)^ (2428691599 + 313957879 + 1279504693 e'' 



+63843717 e^) + 192 (1 + ef (908960573673 - 155717471796 
-88736969547 - 293676299040 - 195313674237 - 26635698156 e^" 
-3467992016^^) y-384(l + ef (336063804453 - 53956775638 
-33318942522 - 92857670352 - 41764459155 - 2765710514 e^°) 
+16 (1 + ef (3418907055555 - 580720618635 - 168432860626 
-606890963686 - 176495184865 - 3768291999 e^°) 

-32 (1 + ef (510454645597 - 92175635794 + 26432814256 - 28250211070 

-5713846269 e^) y"* + 4 (1 + e)*^ (1107402703901 - 174239346926 
+100957560852 + 3707280110 - 899162673 e«) y'^ 

-8 (1 + e)^ (143625217397 - 160328200106^ + 4238287541 + 275190560 e^) y 
+ (1 + e)* (220627324753 - 14884378223 - 1210713997 e"* + 14138955 e^) 
-8 (1 + ef (2922108518 - 46504603 - 2407656 e^) / 
+3 (1 + e)^° (241579935 + 6314675 - 149426 e^) 

+4(l + e)^^ (8608805 - 48992 e^) yi" + 2(l + e)^^ (1242083 - 16320 e^) 
+184320 (1 + e)" _^ 53^20 (1 + ef^ y^^ 
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and 

f2{y, e) = 3072 (3 - e) (3 + e) (3 + e^) (7286074797 - 3299041125 + 792940362 

- 1366777698 e^- 369698151 e^- 5932745 e^°) -384 (1 + e) (2989180413711 

-583867932642 - 131661872359 - 419423580924 - 194293515951 
-3390301442 6^"+ 1353430119 6^2^ y + 64(l + e)^ (14825178681327 
-2675442646782 - 728511901515 - 1837874368340 - 591999524567 
-1856757710 e^° + 841581651 e^2) - 32 (1 + e)^ (14292163934541 
-2666166422089 - 522582885086 - 1347373382962 - 307066297439 
-1675056789 e^°) + 16 (1 + e)^ (9557748374919 - 1917809903861 
-24258045506 - 511875047746 - 86779453317 - 462078345 e^°) 
-8 (1 + e)^ (5390797838491 - 990602472036 + 161182699002 
-89978894004 - 11363685245 e^) + 4(1 + ef (2857676457065 
-351292910556 + 79840371470 - 2670080940 - 463345647 e^) 
-2 (1 + ef (1249768416047 - 79903103833 + 12179840133 
+482157413 e^) / + (1 + e)^ (363565648057 - 10040939153 - 318841465 
+14611473 e*^') y^-2{l + ef (13862653487 - 100645509 - 11015842 e'^) y^ 
+ (l+e)^" (518128485+ 16345427 62- 421398 6"*) 

+16 (1 + 6)" (1220639 - 13448 e^) + 2 (1 + 6)^^ (689123 - 18880 e^) y^^ 
+153600 (1 + 6)*^ + 5120(1 + ef^ y" 

The angular momentum lost is similarly given by 



^ 24249225 (1 + e)^ (fp - 2 M)2 ((1 + e)rp - 2(1 - e)Mf 



(1 + 6) J - 2(3 - e) E ( e) 



(1 + 6) 46 



^ \/((l + eK-2(3-e)M) ^'(m''') 



^(l + e)^_2(3-e) V V + " - 
where 

yi(y,6)= 169728 (1 - 6) (1 + 6)^ (279297 + 2198976^ + 106299 6^ + 96116^) 

-384 (1 + ef (192524061 - 13847615 - 36165965 6* - 20710173 - 588532 6^) y 

+192 (1 + 6)^ (235976417 + 13109547 - 3369705 6^ - 3292707 e^) y^ 

-16(1 + 6)^ (813592799 +112906199 6^ + 53843933 6^ + 602061 6*^) y^ 

+16 (1 + ef (87491089 + 7247482 6^ + 4608349 e^) y^ + 8 (1 + 6)^ (9580616 

+6179243 6^ - 92047 6^) y^ - 4 (1 + 6)^ (3760123 + 272087 6^) y^ 

-(1 + 6)^ (1168355- 35347 6^) y^ - 71792 (1 + 6)® y^ - 4120 (1 + 6)^° y^ 



(A4) 
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and 

£,2(2/, e) = 339456 (3 - e) (3 + e) (93099 - 10213 - 18155 - 10551 - 420 e*) 

-1536 (1 + e) (319648410 - 35712133 - 33099777 - 11272311 + 457187 e^) y 
+ 128 (1 + ef (2706209781 - 45415294 - 103634296 - 34056010 - 130293 e^) 
-32 (1 + e)^ (3895435659 + 212168215 + 4641265 - 15197651 e^) 
+ 16(1 + ef (1396737473 + 123722895 + 27602127 - 465119 e*^) y"* 
-16(l + e)^ (78148621 + 3035912 6^ + 3130827 6"*) y^ 

-16 (1 + ef (8005570 + 1485159 - 47943 e'^) + 2 (1 + ef (4015181 + 601959 e^) 
+ (1 + 6)^ (737603 - 39467 e^) / + 47072 (1 + 6)^ y^ + 4120 (1 + 6)^" y^" 

The limit rp — > cx) corresponds to the argument of the elliptic integrals approaching zero. Using series expansions of 
the elliptic integrals about k = Q [2^], we find 

M^^^_64.^^^ 73. 37 
m 5 (i + e)^ 



24 ^ 96 / \m) 



64 n 



-^('1+^6^+^6^+M f^r^+o 



(1+6) 

m 5 (i + ef V 8 y Vm^ 

192 TT 1 / 35 , 1 , 
1 + :rT + - 6^ 



(1 + e)^ V 24 4 



m Mm ■ <«) 



The leading order terms agree with the Keplerian results ((22)l - ((23l) as expected. In the limit — > rjjco = 
2(3 + e)/(l + e), the argument of the elliptic integrals approaches 1. The elliptic integrals diverge logarithmically in 
this limit, and we may expand them as in equations ([38| and p9|) . 

On substitution of these expansions into (IA3p - (|A4p . we find the asymptotic form of AE and ALz to be 

^AX.p.lnf^-^i^U,..+0ff^-^(^Vnf^-^^)) (A6) 
m \M 1 + 6 y ^ \\M 1 + e J \M l + e JJ ^ ' 

As before, X refers to either E or L^/M . The coefficients px and qx are functions of eccentricity 
5 Ve (3 + e) 

= 4 ^^ (126493657290+ 548139181590 6+ 1030019780790 6^ 

+1139255611065 + 838466930873 + 401719467929 + 98700067049 
+6236043751 6^ + 2856045401 - 177251547 6^ - 1203124043 e^° 



+316812556 6" + 109455696 e^^ - 88995328 e^^) / (^1673196525 (1 + e) ^ (3 + e 
4 (1 + ef (ln(64) - ln(l + 6) + ln{e)) 



5Vi(3 + 6) 



3 



(A8) 



5(3 + e)2 



16V2(l + 6) 



2 



6 



gi, = — '- J ^(174594420+ 523783260 6 

24249225 (3 + 6)^ V + + e) 

+557732175 e^ + 241337525 + 44249062 e'^ + 11244922 e^ - 2993241 6^ - 1809123 6^' 

.o<.„„n. 8 qx 4849845 (61n(2)-ln(l + 6)+ln(6))\ 
+1328784 6^ - 172744 6^) ^ r, — (AlO) 
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2. Fitting functions 

We can use the exact expressions (|A3[) ~ (|A4p to derive fitting functions to approximate our results. Following the 
same argument used in the parabolic case, a functional form like (j29p should capture the main features of the problem, 
but the coefficients are now functions of eccentricity, and we replace the parabolic value of the UCO - 4 Af - with the 
value appropriate to other eccentricities. The general ansatz is 



m 




(1 + e)r2 ) 



BP 



2(3 + e)A/\^^"^ (l + e)M 



" V (1 + e)rp ) (1 + e)rp - 2(3 + e)M 

Af^((l + e)rp-2(3 + e)M) ^ ^ / Af ((1 + e)rp - 2(3 + e)M) 

(1 + ejTp ^ n=0 ^ \ / p 

Afi+^((1 + e)rp - 2(3 + e)M) ^ Af ((1 + e)rp - 2(3 + e)A/) V 

+ ] — nr^i ^ I fi+pV2 ) (All) 

Successive terms in the fit are given by matching consecutive orders in an expansion about Vp = 2(3 + e)/(l + e) and 
as Tp ^ oo in the way described in section IIIID II for the parabolic case. To illustrate, the lowest order {N — 0) 
expansion coefficients may be determined from the AE and AL^ expansions (|A5|) - (|A10p as follows 



^?(e) = -Pxie), B^{e) = - exp 



2 " V^x(e) 



.E,^, 64^ 1 A , 73 , , 37 A /:;^i7:^2(3 + e) 



- --^1^ + ^^ +96^ j-^o(^V2^o^(^)V(l + e) 
_ 64^ 1 (. , 7 ,\ .L...^./:;^!— ;/2(3 + e) 



Co-(e) = ^1+ ^e^j - A-(e)V2i?o-W j (A12) 

Our main focus is on orbits that are nearly parabolic, with e ~ 1. We therefore expand these expressions about e = 1 
to obtain 

B^{e) = 0.752091- 0.0949439(1 -e)+ 0.0918458(1 -e)2+0((l-e)3) 
C^(e) = -4.63464+ 1.63944(1 -e)+ 0.327505(1 -e)2 + 0((l-e)3) 

¥+^<'-)- 

S^-(e) = 1.31899- 0.126207(1 -e) + 0.392812(1 -e)2 + 0((l-e)3) 

Co ^ (e) = -4.1491 + 1.71517(1 - e) - 0.128645(1 - + O ((1 - e)^) (A13) 

The expansion of the B^^s and C(^'s may also be written down precisely. However, the expressions are extremely 
complicated, which is why the numerical values of these coefficients have been quoted instead. Higher order fitting 
functions may be obtained by matching more terms in the expansions of AE and ALz, as described earlier. For 
completeness, we quote here the remaining coefficients of the N = 2 fit, once again expanded to quadratic order about 
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the parabohc case 

Af{e) = -0.282843(1 - e) + 0.0353553(1 - e)^ + O ((1 - e)^) 

Bf(e) = -103.215 + 39.6287(1 -e)+ 38.3325(1 -e)2 + 0((l-e)3) 

Cf(e) = 69.1683- 0.682028(1 -e)- 28.7945(1 -e)2 + 0((l-e)3) 

(e) = -1.20797 - 2.31872(1 - e) - 2.15134(1 - e)^ + O ((1 - e)^) 

Bf(e) = 727.515+ 1570.89(1 -e) + 1139.13(1 -e)2 + 0((l-e)3) 

Cf(e) = -439.378- 1223.38(1 -e) - 862.812(1 -e)2 + 0((l -e)3) 

A^'{e) = -0.565685(1 -e) + 0.494975(1 -e)2 + 0((l-e)^) 

B^'ie) = -53.4491 + 4.38709(1 -e) + 0.469838(1 -e)2 + C>((l-e)3) 

Cf-(e) = 25.4129+ 16.7694(1 -e)- 7.06419(1 -e)2 + 0((l-e)'') 

A^'ie) = 3.9598(1 -e)- 4.80833(1 -e)2 + 0((l-e)3) 

B^'{e) = 29.7857+ 167.281(1 -e) + 66.0607(1 -e)^+0((l-e)^) 

C2'{e) = 15.1726- 131.512(1 -e)- 26.8611(1 -e)^+0((l-e)^) (A14) 



APPENDIX B: ERRATUM — PHYS. REV. D 74 109901 (2006) 



Erratum: Semi-relativistic approximation to gravitational radiation from encounters 

with non-spinning black holes 
[Phys. Rev. D 72, 084009 (2005)] 

Jonathan R Gair, Daniel J Kennefick, Shane L Larson 

In Section III.E on hyperbolic orbits wc stated "The energy and angular momentum lost to gravitational waves by a 
particle on a hyperbolic orbit are given by the same Eqs. (A3)-(A4) that apply to bound orbits, just by inserting e > 1 
consistent with the definitions (12) and (13)". In fact, the generalisation to hyperbolic orbits is more complicated, 
since the integrals are now incomplete, and I2 defined by equation (Al) changes since a boundary term no longer 
cancels. The correct generalisation to e > 1 requires modifying expressions (A3) and (A4) as follows 

• Replace K(A;) by F{^,k), an incomplete elliptic integral of the first kind, where 



.2 _ 4eM ^ = .ir,-M (l + e)((l + e)rp-2(3-e)M) \ 

(l + e)rp-2(3-e)M' ^ []] 2e((l + e)rp - 4M) J 

• Replace E(A;) by E(^,A;) - EH(rp,e), where E((^,fc) is the incomplete elliptic integral of the second kind and 



EH(rp,e) =2M, 



e2-l 



((1 + e)rp - 2(3 - e)M) ((1 + e)rp - 4M) 

• Include an extra term, (1 + e) Ve^ — 1 f3{rp/M, e)/-\/(l + e)rp/M — 4, inside the square bracket of equation 
(A3), and a term, (1 + e) Ve^'^^ g3{rp/M, e)/ ^/JT+'e)r^J^^^^, inside the square bracket of equation (A4). 
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The new functions are given by 

f^{y, e) = 5120(e + l^^y^^ + 186880(e - l)(e + 1)^^?/^^ - 2(e + 1)^^ (l5040e^ - 1294563) y^^ 

-(e + 1)1° (897286^ - 36202743) - 2(e + 1)^ (229899e^ - 106856626^ - 375412235) 
+ {e + if (I9031885e'^ - 2703503366e2 - 25681442087) ?/ + 4(e + 1)^ (3474054e^ - 1305675539e'' 
+3107359416e2 + 59447608277) - 4(e + 1)*^ (865886773e*^' - 22415704847e'' - 304076054096^ 
+296614649595) y*^ + 16(e + 1)^ (534461702e^ + 115330338976*^ - 22057127975e'' - 766037700536^ 
+272389989629)1/^ - 16(e + l)'^ (l7178999909e^ + 162890646772e^ + 592951251586"* 
-2872238849886^ + 971850401469) y^ + 64(6 + 1)^ (6036828186*° + 454062729116^ 
+2432085070406^ + 155132151938e'' - 157983734058e^ + 799957785207) y^ 
-64(e + if (80581457316*° + 210515698689e^ + 7069861723826*^ + 411976993282e'' 
-2636689345136^ + 1861696008653) y"^ + 256(e + 1) (355964706*^ + 8535177147e*° 
+1113176816956^ + 249751373150e^ + 118284328500e'' - 882707366176^ + 620968567495) y 
-2048 (177982356*2 + 14458813446*° + 111312311236^ + 174666676406^ + 64648850736^ 
-78385378486^ + 44213644993) 

93{y, e) = -4120(e + l)''^^ - 73852(e + 1)**/ + (e + 1)^ (332876^ - 1214551) + 4(e + 1)*^ (9297036^ 
+3916957) / + 16(e + 1)^ (2579556^^ + 21255556^ - 10535542) - 16(e + l)'^ (62098676^ 
+364220106^ - 127264229) y^ - 64(e + 1)^ (5324956° + 7491076"^ - 346052916^ + 203477737) y^ 
+128(6 + if (75918716^ + 40714475e^ - 24139707e^ + 328231985) _ 256(6 + 1) (278460e* 
+227377616^ + 79099653e* - 77682856^ + 262081211) y + 226304 {e^ + 3) (l260e*^ + 40039e^ 
-233786^ + 62719) 

The results shown in Figure 8 were computed by numerical integration of the fluxes, and are therefore correct inspite 
of this error. This erratum has been published as Phys. Rev. D 74 109901. 



